2020 Environmental Performance Index (EPI)

“The 2020 Environmental Performance Index (EPI) provides a data-driven summary of the state of sustainability around the world. Using 32 performance indicators across 11 issue categories, the EPI ranks 180 countries on environmental health and ecosystem vitality. These indicators provide a gauge at a national scale of how close countries are to established environmental policy targets. The EPI offers a scorecard that highlights leaders and laggards in environmental performance and provides practical guidance for countries that aspire to move toward a sustainable future.” – 2020 Environmental Performance Index

Load libraries

library("reshape2")
library("ggplot2")
theme_set(theme_bw(base_size = 12) +
            theme(rect = element_rect(fill = "transparent")))
library("leaflet")
library("rgdal")
library("htmlwidgets")
library("lavaan")
library("lavaanPlot")

Read data

Clear workspace, set working directory

rm(list = ls())
setwd("~/epi")

Load helper functions

These helper functions retrieve, read and process data.

source("UtilityFunctions.R")

Import EPI data

Run the helper function importEPI() and release the contents of the resulting list into the global environment.

epi.list <- importEPI()
## Warning in rm(epi, envir = .GlobalEnv): object 'epi' not found
## 'data.frame':    8234 obs. of  7 variables:
##  $ code         : int  4 24 8 784 32 51 28 36 40 31 ...
##  $ iso          : chr  "AFG" "AGO" "ALB" "ARE" ...
##  $ country      : Factor w/ 179 levels "Afghanistan",..: 1 4 2 170 6 7 5 8 9 10 ...
##  $ region       : Factor w/ 8 levels "Asia-Pacific",..: 7 8 2 5 6 3 6 4 4 3 ...
##  $ X2019        : num  1.93e+10 8.88e+10 1.53e+10 4.21e+11 4.45e+11 ...
##  $ EPI.new      : Factor w/ 46 levels "AGR","AIR","APE",..: 12 12 12 12 12 12 12 12 12 12 ...
##  $ EPI.new.value: num  25.5 29.7 49 55.6 52.2 52.3 48.5 74.9 79.6 46.5 ...

list2env(epi.list, .GlobalEnv)
## <environment: R_GlobalEnv>
rm(epi.list)

Download world borders

Download the shape file, if it is absent. Note, this overwrites existing files!

if(!file.exists("2020/TM_WORLD_BORDERS-0.3.shp")){
  print("File does not exist, downloading it now:")
  download.file("thematicmapping.org/downloads/TM_WORLD_BORDERS-0.3.zip",
                destfile="world_shape_file.zip")
  system("unzip -u world_shape_file.zip;
         rm world_shape_file.zip Readme.txt")
  }

Then, read this shape file into a SpatialPolygonsDataFrame and inspect it.

world.spdf <- readOGR(dsn = "2020", layer = "TM_WORLD_BORDERS-0.3",
                      verbose = TRUE)

## OGR data source with driver: ESRI Shapefile 
## Source: "/home/r/epi/2020", layer: "TM_WORLD_BORDERS-0.3"
## with 246 features
## It has 11 fields
## Integer64 fields read as strings:  POP2005

summary(world.spdf)

## Object of class SpatialPolygonsDataFrame
## Coordinates:
##    min      max
## x -180 180.0000
## y  -90  83.6236
## Is projected: FALSE 
## proj4string : [+proj=longlat +datum=WGS84 +no_defs]
## Data attributes:
##      FIPS               ISO2               ISO3                 UN       
##  Length:246         Length:246         Length:246         Min.   :  4.0  
##  Class :character   Class :character   Class :character   1st Qu.:215.0  
##  Mode  :character   Mode  :character   Mode  :character   Median :429.0  
##                                                           Mean   :431.8  
##                                                           3rd Qu.:650.5  
##                                                           Max.   :894.0  
##      NAME                AREA             POP2005              REGION      
##  Length:246         Min.   :      0.0   Length:246         Min.   :  0.00  
##  Class :character   1st Qu.:     44.5   Class :character   1st Qu.:  2.00  
##  Mode  :character   Median :   5515.5   Mode  :character   Median : 19.00  
##                     Mean   :  52696.1                      Mean   : 65.43  
##                     3rd Qu.:  34708.8                      3rd Qu.:142.00  
##                     Max.   :1638094.0                      Max.   :150.00  
##    SUBREGION           LON               LAT          
##  Min.   :  0.00   Min.   :-178.13   Min.   :-80.4460  
##  1st Qu.: 14.00   1st Qu.: -50.16   1st Qu.: -0.3025  
##  Median : 30.00   Median :  17.66   Median : 16.5110  
##  Mean   : 54.84   Mean   :  13.29   Mean   : 16.4075  
##  3rd Qu.: 61.00   3rd Qu.:  50.01   3rd Qu.: 39.1067  
##  Max.   :155.00   Max.   : 179.22   Max.   : 78.8300

Collate countries

Unfortunately, as soon as it comes to map the world, geopolitics come to the fore. This is reflected in the mismatches between country lists due to the different recognition statuses of these countries. Luckily, an ISO standard exists that allocates a unique 3-letter code to every country.

For one, there are distinctly more countries in the shp file:

length(world.spdf$NAME) 
## [1] 246
length(epi$country)
## [1] 179

These mismatches stem mainly from different naming schemes…

sort(setdiff(world.spdf$NAME, epi$country))

##  [1] "Åland Islands"                            
##  [2] "American Samoa"                           
##  [3] "Andorra"                                  
##  [4] "Anguilla"                                 
##  [5] "Antarctica"                               
##  [6] "Aruba"                                    
##  [7] "Bermuda"                                  
##  [8] "Bouvet Island"                            
##  [9] "British Indian Ocean Territory"           
## [10] "British Virgin Islands"                   
## [11] "Burma"                                    
## [12] "Cape Verde"                               
## [13] "Cayman Islands"                           
## [14] "Christmas Island"                         
## [15] "Cocos (Keeling) Islands"                  
## [16] "Congo"                                    
## [17] "Cook Islands"                             
## [18] "Democratic Republic of the Congo"         
## [19] "Falkland Islands (Malvinas)"              
## [20] "Faroe Islands"                            
## [21] "French Guiana"                            
## [22] "French Polynesia"                         
## [23] "French Southern and Antarctic Lands"      
## [24] "Gibraltar"                                
## [25] "Greenland"                                
## [26] "Guadeloupe"                               
## [27] "Guam"                                     
## [28] "Guernsey"                                 
## [29] "Heard Island and McDonald Islands"        
## [30] "Holy See (Vatican City)"                  
## [31] "Hong Kong"                                
## [32] "Iran (Islamic Republic of)"               
## [33] "Isle of Man"                              
## [34] "Jersey"                                   
## [35] "Korea, Democratic People's Republic of"   
## [36] "Korea, Republic of"                       
## [37] "Lao People's Democratic Republic"         
## [38] "Libyan Arab Jamahiriya"                   
## [39] "Liechtenstein"                            
## [40] "Macau"                                    
## [41] "Martinique"                               
## [42] "Mayotte"                                  
## [43] "Micronesia, Federated States of"          
## [44] "Monaco"                                   
## [45] "Montserrat"                               
## [46] "Nauru"                                    
## [47] "Netherlands Antilles"                     
## [48] "New Caledonia"                            
## [49] "Niue"                                     
## [50] "Norfolk Island"                           
## [51] "Northern Mariana Islands"                 
## [52] "Palau"                                    
## [53] "Palestine"                                
## [54] "Pitcairn Islands"                         
## [55] "Puerto Rico"                              
## [56] "Republic of Moldova"                      
## [57] "Reunion"                                  
## [58] "Saint Barthelemy"                         
## [59] "Saint Helena"                             
## [60] "Saint Kitts and Nevis"                    
## [61] "Saint Martin"                             
## [62] "Saint Pierre and Miquelon"                
## [63] "San Marino"                               
## [64] "Somalia"                                  
## [65] "South Georgia South Sandwich Islands"     
## [66] "Svalbard"                                 
## [67] "Swaziland"                                
## [68] "Syrian Arab Republic"                     
## [69] "Taiwan"                                   
## [70] "The former Yugoslav Republic of Macedonia"
## [71] "Tokelau"                                  
## [72] "Turks and Caicos Islands"                 
## [73] "Tuvalu"                                   
## [74] "United Republic of Tanzania"              
## [75] "United States"                            
## [76] "United States Minor Outlying Islands"     
## [77] "United States Virgin Islands"             
## [78] "Wallis and Futuna Islands"                
## [79] "Western Sahara"                           
## [80] "Yemen"

… but also from the fact that the EPI is missing for some countries. Countries without an EPI are mostly miniature states (e.g. Monaco, Liechtenstein or Andorra), states currently in a(n) (armed) conflict (Syria, Palestine or Taiwan) or several islands states. Some of these islands actually belong to a nation; e.g., Greenland to Denkmark and Svalbard to Norway (for now, keep it as is although it would be reasonable to merge these islands with their mainlands).

sort(setdiff(epi$country, world.spdf$NAME))

##  [1] "Cabo Verde"               "Dem. Rep. Congo"         
##  [3] "Eswatini"                 "Iran"                    
##  [5] "Laos"                     "Micronesia"              
##  [7] "Moldova"                  "Myanmar"                 
##  [9] "North Macedonia"          "Republic of Congo"       
## [11] "South Korea"              "Tanzania"                
## [13] "United States of America"

Merge datasets

Now, merge the epi data.frame and the world.spdf SpatialPolygonsDataFrame by the country ISO-3 so that the latter holds all relevant information to render the map.

world.spdf <- merge(world.spdf,
                    epi[, (colnames(epi) %in% c("iso", "EPI.new", "region",
                                                "X2019"))], 
                    by.x = "ISO3", by.y = "iso")
names(world.spdf)

##  [1] "ISO3"      "FIPS"      "ISO2"      "UN"        "NAME"      "AREA"     
##  [7] "POP2005"   "REGION"    "SUBREGION" "LON"       "LAT"       "EPI.new"  
## [13] "region"    "X2019"

Encode population numbers

For readability, population numbers are encoded as millions (M) and for meaningfulness, as numbers instead of characters.

world.spdf@data$POP2005[which(world.spdf@data$POP2005 == 0)] <- NA
world.spdf@data$POP2005 <- as.integer(world.spdf@data$POP2005) / 1000000 %>%
  round(2)

Choropleth map

Now, let’s render an interactive map with R leaflet to display the EPI on the globe.

Color palette

First, create a gradient color scheme for the EPI:

epi.cols <- colorNumeric(palette = "viridis", domain = world.spdf@data$EPI.new,
                         na.color = "gray", reverse = TRUE)

An then a discrete color scheme for the regions:

region.cols <- colorFactor(palette = "plasma", domain = world.spdf@data$region,
                           na.color = "gray")

Tooltips

These tooltips emerge when hovering over the html-ified map showing some data of the country under the mouse pointer.

tooltips <- paste("Country:", world.spdf@data$NAME, "<br/>",
                  "Region:", world.spdf@data$region, "<br/r>",
                  "GDP:", round(world.spdf@data$X2019 / 1000000000, digits = 1),
                  "B US$", "<br/r>",
                  # "Area:", world.spdf@data$AREA, "<br/>",
                  "Population:", round(world.spdf@data$POP2005, 2), "M",
                  "<br/>",
                  "EPI:", world.spdf@data$EPI.new) %>%
  lapply(htmltools::HTML)

Render map

Finally, use the previous steps to generate a leaflet choropleth map.

epi.map <- leaflet(world.spdf) %>%
  addTiles() %>%
  setView(lat = 20, lng = 0, zoom = 1.5) %>%
  addPolygons(fillColor = ~ epi.cols(EPI.new), stroke = TRUE, group = "EPI",
              fillOpacity = 0.9, color = "white", weight = 0.5,
              label = tooltips,
              labelOptions = labelOptions(style = list("font-weight" = "normal",
                                                       padding = "3px 8px"),
                                          textsize = "13px",
                                          direction = "auto")) %>%
  addPolygons(fillColor = ~ region.cols(region), stroke = TRUE,
              group = "Region", fillOpacity = 0.5, color = "white",
              weight = 0.5, label = tooltips,
              labelOptions = labelOptions(style = list("font-weight" = "normal",
                                                       padding = "3px 8px"),
                                          textsize = "13px",
                                          direction = "auto")) %>%
  addLegend(pal = epi.cols, values = ~ EPI.new, opacity = 0.9,
            title = "EPI", position = "topleft", group = "EPI") %>%
  addLegend(pal = region.cols, values = ~ region, opacity = 0.5,
            title = "", position = "topleft", group = "Region") %>%
  addLayersControl(overlayGroups = c("EPI", "Region"),
                   options = layersControlOptions(collapsed = FALSE)) %>%
  hideGroup("Region")
epi.map

Export map to html

The map can easily be exported to html.

saveWidget(epi.map, file = "index.html")

Analysis

Color scale for plotting

For consistency, use the same colors for the regions as before in the choropleth map.

cols.region <- c("#0D0887", "#5402A3", "#8B0AA5", "#B93289", "#DB5C68",
                 "#F48849", "#FEBC2A", "#F0F921")

National comparison

To obtain an overview of the worldwide distribution of the EPI, the 32 performance indicators are displayed by country in boxplots sorted by their median and colored by region.

ggplot(epi.long[epi.long$Type == "Indicator", ],
       aes(x = reorder(country, EPI.new.value, FUN = median, na.rm = TRUE),
           y = EPI.new.value, fill = region)) +
  geom_boxplot(outlier.shape = 21, outlier.size = 1, alpha = 0.5) +
  scale_fill_manual(values = cols.region) +
  theme(legend.position = "top", legend.title = element_blank(),
        axis.text.x = element_text(angle = 45, hjust = 1),
        axis.text = element_text(size = 5)) +
  ylab("2020 performance indicators") +
  xlab("")
## Warning: Removed 446 rows containing non-finite values (stat_boxplot).

# ggsave("EPI_nations.pdf", width = 20, height = 8.27, bg = "transparent")

An ordering according to region is apparent. So, let’s have a closer look onto these regions.

Regional comparison

To highlight the regional differences, the EPIs are pooled by region and displayed in boxplots sorted by their median and colored by region. That way, the strinking gap between the Global West to the remaining regions becomes apparent.

ggplot(epi.long[epi.long$Type == "EPI", ],
       aes(x = reorder(region, EPI.new.value, FUN = median, na.rm = TRUE),
           y = EPI.new.value, fill = region)) +
  geom_boxplot(outlier.shape = 21, outlier.size = 1, alpha = 0.5) +
  scale_fill_manual(values = cols.region) +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1)) +
  ylim(0, 100) +
  ylab("2020 EPI") +
  xlab("")

# ggsave("EPI_regions.pdf", width = 11.29, height = 6.27, bg = "transparent")

EPI vs GDP

There seems to be a distinct relation between the wealth of a country and its EPI. Plotting the gross domestic product against the EPI shows that this actually true. However, this is deceptive as the GDP is on a logarithmic scale. The correlation is actually only 0.20, but increases to 0.53 with the log scaled GDP.

ggplot(epi.long[epi.long$Type == "EPI", ],
       aes(x = log(X2019 / 1000000), y = EPI.new.value, fill = region)) +
  geom_smooth(aes(group = 1), color = "gray", lty = 2, method = "lm",
              formula = y ~ x, se = FALSE, show.legend = FALSE) +
  geom_point(alpha = 0.5, pch = 21) +
  scale_fill_manual(values = cols.region) +
  theme(legend.position = "top", legend.title = element_blank()) +
  ylim(0, 100) +
  ylab("2020 EPI") +
  xlab(expression(paste(log[e], "(2019 GDP)", " [US$]")))
## Warning: Removed 7 rows containing non-finite values (stat_smooth).
## Warning: Removed 7 rows containing missing values (geom_point).

# ggsave("EPI_GDP.pdf", width = 11.29, height = 6.27, bg = "transparent")

Structural equation model

Now, let’s use a structural equation model to reveal the driving forces behind the EPI.

Merge datasets

epi <- merge(epi,
             world.spdf[, (names(world.spdf) %in% c("ISO3", "AREA",
                                                    "POP2005", "LAT"))],
             by.x = "iso", by.y = "ISO3")

Recode latitude

For a straight-forward comparison, the latitude is reset to 0 at the south pole and 180 at the north pole by adding 90.

epi$LAT <- epi$LAT + 90

Scale the variables

For a SEM, the variables need to have a similar variance. Thus, the \(z\)-scores are calculated to bring the variables on a mean of zero and a standard deviation of one.

epi[, c(4, 51:54)] <- apply(epi[, c(4, 51:54)], 2, scale)
summary(epi[, c(4, 51:54)])

##     EPI.new            X2019              AREA             POP2005       
##  Min.   :-1.5371   Min.   :-0.2418   Min.   :-0.36828   Min.   :-0.2545  
##  1st Qu.:-0.7907   1st Qu.:-0.2359   1st Qu.:-0.35438   1st Qu.:-0.2447  
##  Median :-0.1541   Median :-0.2224   Median :-0.30077   Median :-0.2126  
##  Mean   : 0.0000   Mean   : 0.0000   Mean   : 0.00000   Mean   : 0.0000  
##  3rd Qu.: 0.5052   3rd Qu.:-0.1181   3rd Qu.:-0.09043   3rd Qu.:-0.1067  
##  Max.   : 2.3342   Max.   :10.2255   Max.   : 8.21141   Max.   : 9.5908  
##                    NA's   :7                                             
##       LAT          
##  Min.   :-2.52781  
##  1st Qu.:-0.61659  
##  Median :-0.06604  
##  Mean   : 0.00000  
##  3rd Qu.: 0.87677  
##  Max.   : 1.89993  
## 

Specify the model

The assumption is that the EPI is positivly driven by the GDP, the population size and the latitude with following reasoning:

• the GDP as a measure of economic development is kinda prerequisite for environmental development
• larger population size is expected to result in a higher GDP due to higher humanpower
• the latitude represents world’s north-south divide

The GDP, in turn, is influenced by population size and latitude and the population size by a country’s area.

epi.gdp <- 'EPI.new ~ X2019 + POP2005 + LAT
  X2019 ~ LAT + POP2005
  POP2005 ~ AREA'

Run the model and inspect the results

Running the model reveals a reasonable fit for the robust estimator.

fit.epi.gdp <- sem(epi.gdp, data = epi, estimator = "MLM")
summary(fit.epi.gdp, rsq = TRUE)

## lavaan 0.6-7 ended normally after 13 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of free parameters                          9
##                                                       
##                                                   Used       Total
##   Number of observations                           172         179
##                                                                   
## Model Test User Model:
##                                               Standard      Robust
##   Test Statistic                                26.066       5.408
##   Degrees of freedom                                 3           3
##   P-value (Chi-square)                           0.000       0.144
##   Scaling correction factor                                  4.820
##        Satorra-Bentler correction                                 
## 
## Parameter Estimates:
## 
##   Standard errors                           Robust.sem
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   EPI.new ~                                           
##     X2019             0.292    0.127    2.311    0.021
##     POP2005          -0.271    0.073   -3.727    0.000
##     LAT               0.509    0.079    6.457    0.000
##   X2019 ~                                             
##     LAT               0.120    0.047    2.525    0.012
##     POP2005           0.582    0.203    2.872    0.004
##   POP2005 ~                                           
##     AREA              0.459    0.250    1.838    0.066
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .EPI.new           0.657    0.074    8.897    0.000
##    .X2019             0.622    0.466    1.335    0.182
##    .POP2005           0.815    0.476    1.714    0.087
## 
## R-Square:
##                    Estimate
##     EPI.new           0.357
##     X2019             0.372
##     POP2005           0.211

Inspect modification indices

The modification indices suggest–among others–that also area influences the GDP. This makes sense as a larger country also has the potential to posess more resources. Thus, this path is added to the model what further increases the model fit.

modificationindices(fit.epi.gdp, minimum.value = 3)
##        lhs op     rhs     mi    epc sepc.lv sepc.all sepc.nox
## 15   X2019 ~~ POP2005 24.039 -0.581  -0.581   -0.816   -0.816
## 18   X2019  ~    AREA 24.039  0.327   0.327    0.334    0.329
## 20 POP2005  ~   X2019 20.315 -0.814  -0.814   -0.797   -0.797
## 27    AREA  ~   X2019 21.889  0.494   0.494    0.484    0.484
fit.epi.gdp.up <- update(fit.epi.gdp, add = "X2019 ~ AREA")
summary(fit.epi.gdp.up, rsq = TRUE)
## lavaan 0.6-7 ended normally after 16 iterations
## 
##   Estimator                                         ML
##   Optimization method                           NLMINB
##   Number of free parameters                         10
##                                                       
##                                                   Used       Total
##   Number of observations                           172         179
##                                                                   
## Model Test User Model:
##                                               Standard      Robust
##   Test Statistic                                 0.211       0.204
##   Degrees of freedom                                 2           2
##   P-value (Chi-square)                           0.900       0.903
##   Scaling correction factor                                  1.034
##        Satorra-Bentler correction                                 
## 
## Parameter Estimates:
## 
##   Standard errors                           Robust.sem
##   Information                                 Expected
##   Information saturated (h1) model          Structured
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)
##   EPI.new ~                                           
##     X2019             0.292    0.127    2.310    0.021
##     POP2005          -0.271    0.073   -3.725    0.000
##     LAT               0.509    0.079    6.442    0.000
##   X2019 ~                                             
##     LAT               0.109    0.042    2.615    0.009
##     POP2005           0.432    0.175    2.471    0.013
##   POP2005 ~                                           
##     AREA              0.459    0.250    1.838    0.066
##   X2019 ~                                             
##     AREA              0.327    0.253    1.290    0.197
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)
##    .EPI.new           0.657    0.074    8.897    0.000
##    .X2019             0.535    0.360    1.486    0.137
##    .POP2005           0.815    0.476    1.714    0.087
## 
## R-Square:
##                    Estimate
##     EPI.new           0.358
##     X2019             0.460
##     POP2005           0.211

Visualize the model

Finally, the results are plotted as a path diagram including the path coefficients and their significance levels. This path diagram shows, that based on the available data, the strongest direct influence on the EPI is the latitude followed by the GDP and in a negative way, the population size.

lavaanPlot(model = fit.epi.gdp.up, coefs = TRUE, stars = c("regress"),
           node_options = list(shape = "box", color = "gray",
                               fontname = "Helvetica"),
           edge_options = list(color = "black"),
           labels = list(AREA = "Area", POP2005 = "Population", X2019 = "GDP",
                         LAT = "Latitude", EPInew = "EPI"))

Conclusion

Each step added a further insight onto the EPI:

1. The choropleth map showed the distribution of high- and low-valued countries across the globe, whereby a distinct cluster of high-value countries in North America, Europe and Australia was apparent and lowest values were found in East Asia and Africa
2. This finding was supported by the inspection of the EPI by country and region. Especially the regional display revealed a strinking gap between the Global West to the remaining regions.
3. The idea that the GDP drives the EPI seems intuitive: richer countries potentially invest more into a healthy environment. However, the correlation was quite low.
4. The SEM showed that the north-south divide is also reflected in the EPI. Interestingly, the GDP was only the second-most important influence of the EPI. However, this result may likely change when using the purchasing power parity (PPP) instead of the GDP, because large countries such as China or Brazil do well regarding the GDP, but perform comparatively poorly regarding the PPP.

Session info

sessionInfo()

## R version 4.0.4 (2021-02-15)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Debian GNU/Linux 10 (buster)
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.8.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.8.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] lavaanPlot_0.6.0  lavaan_0.6-7      htmlwidgets_1.5.3 rgdal_1.5-23     
## [5] sp_1.4-5          leaflet_2.0.4.1   ggplot2_3.3.3     reshape2_1.4.4   
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.1.0   xfun_0.20          purrr_0.3.4        splines_4.0.4     
##  [5] lattice_0.20-41    colorspace_2.0-0   vctrs_0.3.6        generics_0.1.0    
##  [9] htmltools_0.5.1.1  stats4_4.0.4       viridisLite_0.3.0  mgcv_1.8-33       
## [13] yaml_2.2.1         rlang_0.4.10       pillar_1.4.7       glue_1.4.2        
## [17] withr_2.4.1        DBI_1.1.1          RColorBrewer_1.1-2 lifecycle_0.2.0   
## [21] plyr_1.8.6         stringr_1.4.0      munsell_0.5.0      gtable_0.3.0      
## [25] visNetwork_2.0.9   evaluate_0.14      labeling_0.4.2     knitr_1.31        
## [29] crosstalk_1.1.1    DiagrammeR_1.0.6.1 highr_0.8          Rcpp_1.0.6        
## [33] scales_1.1.1       jsonlite_1.7.2     tmvnsim_1.0-2      farver_2.0.3      
## [37] gridExtra_2.3      mnormt_2.0.2       digest_0.6.27      stringi_1.5.3     
## [41] dplyr_1.0.4        grid_4.0.4         tools_4.0.4        magrittr_2.0.1    
## [45] tibble_3.0.6       crayon_1.4.1       pbivnorm_0.6.0     pkgconfig_2.0.3   
## [49] Matrix_1.3-2       ellipsis_0.3.1     rstudioapi_0.13    assertthat_0.2.1  
## [53] rmarkdown_2.6      viridis_0.5.1      R6_2.5.0           nlme_3.1-152      
## [57] compiler_4.0.4